home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 November / macformat-018.iso / Utility Spectacular / Utilities / Calc / zmod.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-10-10  |  31.7 KB  |  1,338 lines  |  [TEXT/????]
  1. /*
  2.  * Copyright (c) 1992 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Routines to do modulo arithmetic both normally and also using the REDC
  7.  * algorithm given by Peter L. Montgomery in Mathematics of Computation,
  8.  * volume 44, number 170 (April, 1985).  For multiple multiplies using
  9.  * the same large modulus, the REDC algorithm avoids the usual division
  10.  * by the modulus, instead replacing it with two multiplies or else a
  11.  * special algorithm.  When these two multiplies or the special algorithm
  12.  * are faster then the division, then the REDC algorithm is better.
  13.  */
  14.  
  15. #include "xmath.h"
  16.  
  17.  
  18. #define    POWBITS    4        /* bits for power chunks (must divide BASEB) */
  19. #define    POWNUMS    (1<<POWBITS)    /* number of powers needed in table */
  20.  
  21.  
  22. LEN _pow2_ = POW_ALG2;        /* modulo size to use REDC for powers */
  23. LEN _redc2_ = REDC_ALG2;    /* modulo size to use second REDC algorithm */
  24.  
  25. static REDC *powermodredc = NULL;    /* REDC info for raising to power */
  26.  
  27. #if 0
  28. extern void zaddmod proto((ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res));
  29. extern void znegmod proto((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  30.  
  31. /*
  32.  * Multiply two numbers together and then mod the result with a third number.
  33.  * The two numbers to be multiplied can be negative or out of modulo range.
  34.  * The result will be in the range 0 to the modulus - 1.
  35.  */
  36. void
  37. zmulmod(z1, z2, z3, res)
  38.     ZVALUE z1;        /* first number to be multiplied */
  39.     ZVALUE z2;        /* second number to be multiplied */
  40.     ZVALUE z3;        /* number to take mod with */
  41.     ZVALUE *res;        /* result */
  42. {
  43.     ZVALUE tmp;
  44.     FULL prod;
  45.     FULL digit;
  46.     BOOL neg;
  47.  
  48.     if (iszero(z3) || isneg(z3))
  49.         error("Mod of non-positive integer");
  50.     if (iszero(z1) || iszero(z2) || isunit(z3)) {
  51.         *res = _zero_;
  52.         return;
  53.     }
  54.  
  55.     /*
  56.      * If the modulus is a single digit number, then do the result
  57.      * cheaply.  Check especially for a small power of two.
  58.      */
  59.     if (istiny(z3)) {
  60.         neg = (z1.sign != z2.sign);
  61.         digit = z3.v[0];
  62.         if ((digit & -digit) == digit) {    /* NEEDS 2'S COMP */
  63.             prod = ((FULL) z1.v[0]) * ((FULL) z2.v[0]);
  64.             prod &= (digit - 1);
  65.         } else {
  66.             z1.sign = 0;
  67.             z2.sign = 0;
  68.             prod = (FULL) zmodi(z1, (long) digit);
  69.             prod *= (FULL) zmodi(z2, (long) digit);
  70.             prod %= digit;
  71.         }
  72.         if (neg && prod)
  73.             prod = digit - prod;
  74.         itoz((long) prod, res);
  75.         return;
  76.     }
  77.  
  78.     /*
  79.      * The modulus is more than one digit.
  80.      * Actually do the multiply and divide if necessary.
  81.      */
  82.     zmul(z1, z2, &tmp);
  83.     if (ispos(tmp) && ((tmp.len < z3.len) || ((tmp.len == z3.len) &&
  84.         (tmp.v[tmp.len-1] < z2.v[z3.len-1]))))
  85.     {
  86.         *res = tmp;
  87.         return;
  88.     }
  89.     zmod(tmp, z3, res);
  90.     freeh(tmp.v);
  91. }
  92.  
  93.  
  94. /*
  95.  * Square a number and then mod the result with a second number.
  96.  * The number to be squared can be negative or out of modulo range.
  97.  * The result will be in the range 0 to the modulus - 1.
  98.  */
  99. void
  100. zsquaremod(z1, z2, res)
  101.     ZVALUE z1;        /* number to be squared */
  102.     ZVALUE z2;        /* number to take mod with */
  103.     ZVALUE *res;        /* result */
  104. {
  105.     ZVALUE tmp;
  106.     FULL prod;
  107.     FULL digit;
  108.  
  109.     if (iszero(z2) || isneg(z2))
  110.         error("Mod of non-positive integer");
  111.     if (iszero(z1) || isunit(z2)) {
  112.         *res = _zero_;
  113.         return;
  114.     }
  115.  
  116.     /*
  117.      * If the modulus is a single digit number, then do the result
  118.      * cheaply.  Check especially for a small power of two.
  119.      */
  120.     if (istiny(z2)) {
  121.         digit = z2.v[0];
  122.         if ((digit & -digit) == digit) {    /* NEEDS 2'S COMP */
  123.             prod = (FULL) z1.v[0];
  124.             prod = (prod * prod) & (digit - 1);
  125.         } else {
  126.             z1.sign = 0;
  127.             prod = (FULL) zmodi(z1, (long) digit);
  128.             prod = (prod * prod) % digit;
  129.         }
  130.         itoz((long) prod, res);
  131.         return;
  132.     }
  133.  
  134.     /*
  135.      * The modulus is more than one digit.
  136.      * Actually do the square and divide if necessary.
  137.      */
  138.     zsquare(z1, &tmp);
  139.     if ((tmp.len < z2.len) ||
  140.         ((tmp.len == z2.len) && (tmp.v[tmp.len-1] < z2.v[z2.len-1]))) {
  141.             *res = tmp;
  142.             return;
  143.     }
  144.     zmod(tmp, z2, res);
  145.     freeh(tmp.v);
  146. }
  147.  
  148.  
  149. /*
  150.  * Add two numbers together and then mod the result with a third number.
  151.  * The two numbers to be added can be negative or out of modulo range.
  152.  * The result will be in the range 0 to the modulus - 1.
  153.  */
  154. static void
  155. zaddmod(z1, z2, z3, res)
  156.     ZVALUE z1;        /* first number to be added */
  157.     ZVALUE z2;        /* second number to be added */
  158.     ZVALUE z3;        /* number to take mod with */
  159.     ZVALUE *res;        /* result */
  160. {
  161.     ZVALUE tmp;
  162.     FULL sumdigit;
  163.     FULL moddigit;
  164.  
  165.     if (iszero(z3) || isneg(z3))
  166.         error("Mod of non-positive integer");
  167.     if ((iszero(z1) && iszero(z2)) || isunit(z3)) {
  168.         *res = _zero_;
  169.         return;
  170.     }
  171.     if (istwo(z2)) {
  172.         if ((z1.v[0] + z2.v[0]) & 0x1)
  173.             *res = _one_;
  174.         else
  175.             *res = _zero_;
  176.         return;
  177.     }
  178.     zadd(z1, z2, &tmp);
  179.     if (isneg(tmp) || (tmp.len > z3.len)) {
  180.         zmod(tmp, z3, res);
  181.         freeh(tmp.v);
  182.         return;
  183.     }
  184.     sumdigit = tmp.v[tmp.len - 1];
  185.     moddigit = z3.v[z3.len - 1];
  186.     if ((tmp.len < z3.len) || (sumdigit < moddigit)) {
  187.         *res = tmp;
  188.         return;
  189.     }
  190.     if (sumdigit < 2 * moddigit) {
  191.         zsub(tmp, z3, res);
  192.         freeh(tmp.v);
  193.         return;
  194.     }
  195.     zmod(tmp, z2, res);
  196.     freeh(tmp.v);
  197. }
  198.  
  199.  
  200. /*
  201.  * Subtract two numbers together and then mod the result with a third number.
  202.  * The two numbers to be subtract can be negative or out of modulo range.
  203.  * The result will be in the range 0 to the modulus - 1.
  204.  */
  205. void
  206. zsubmod(z1, z2, z3, res)
  207.     ZVALUE z1;        /* number to be subtracted from */
  208.     ZVALUE z2;        /* number to be subtracted */
  209.     ZVALUE z3;        /* number to take mod with */
  210.     ZVALUE *res;        /* result */
  211. {
  212.     if (iszero(z3) || isneg(z3))
  213.         error("Mod of non-positive integer");
  214.     if (iszero(z2)) {
  215.         zmod(z1, z3, res);
  216.         return;
  217.     }
  218.     if (iszero(z1)) {
  219.         znegmod(z2, z3, res);
  220.         return;
  221.     }
  222.     if ((z1.sign == z2.sign) && (z1.len == z2.len) &&
  223.         (z1.v[0] == z2.v[0]) && (zcmp(z1, z2) == 0)) {
  224.             *res = _zero_;
  225.             return;
  226.     }
  227.     z2.sign = !z2.sign;
  228.     zaddmod(z1, z2, z3, res);
  229. }
  230.  
  231.  
  232. /*
  233.  * Calculate the negative of a number modulo another number.
  234.  * The number to be negated can be negative or out of modulo range.
  235.  * The result will be in the range 0 to the modulus - 1.
  236.  */
  237. static void
  238. znegmod(z1, z2, res)
  239.     ZVALUE z1;        /* number to take negative of */
  240.     ZVALUE z2;        /* number to take mod with */
  241.     ZVALUE *res;        /* result */
  242. {
  243.     int sign;
  244.     int cv;
  245.  
  246.     if (iszero(z2) || isneg(z2))
  247.         error("Mod of non-positive integer");
  248.     if (iszero(z1) || isunit(z2)) {
  249.         *res = _zero_;
  250.         return;
  251.     }
  252.     if (istwo(z2)) {
  253.         if (z1.v[0] & 0x1)
  254.             *res = _one_;
  255.         else
  256.             *res = _zero_;
  257.         return;
  258.     }
  259.  
  260.     /*
  261.      * If the absolute value of the number is within the modulo range,
  262.      * then the result is just a copy or a subtraction.  Otherwise go
  263.      * ahead and negate and reduce the result.
  264.      */
  265.     sign = z1.sign;
  266.     z1.sign = 0;
  267.     cv = zrel(z1, z2);
  268.     if (cv == 0) {
  269.         *res = _zero_;
  270.         return;
  271.     }
  272.     if (cv < 0) {
  273.         if (sign)
  274.             zcopy(z1, res);
  275.         else
  276.             zsub(z2, z1, res);
  277.         return;
  278.     }
  279.     z1.sign = !sign;
  280.     zmod(z1, z2, res);
  281. }
  282. #endif
  283.  
  284.  
  285. /*
  286.  * Calculate the number congruent to the given number whose absolute
  287.  * value is minimal.  The number to be reduced can be negative or out of
  288.  * modulo range.  The result will be within the range -int((modulus-1)/2)
  289.  * to int(modulus/2) inclusive.  For example, for modulus 7, numbers are
  290.  * reduced to the range [-3, 3], and for modulus 8, numbers are reduced to
  291.  * the range [-3, 4].
  292.  */
  293. void
  294. zminmod(z1, z2, res)
  295.     ZVALUE z1;        /* number to find minimum congruence of */
  296.     ZVALUE z2;        /* number to take mod with */
  297.     ZVALUE *res;        /* result */
  298. {
  299.     ZVALUE tmp1, tmp2;
  300.     int sign;
  301.     int cv;
  302.  
  303.     if (iszero(z2) || isneg(z2))
  304.         error("Mod of non-positive integer");
  305.     if (iszero(z1) || isunit(z2)) {
  306.         *res = _zero_;
  307.         return;
  308.     }
  309.     if (istwo(z2)) {
  310.         if (isodd(z1))
  311.             *res = _one_;
  312.         else
  313.             *res = _zero_;
  314.         return;
  315.     }
  316.  
  317.     /*
  318.      * Do a quick check to see if the number is very small compared
  319.      * to the modulus.  If so, then the result is obvious.
  320.      */
  321.     if (z1.len < z2.len - 1) {
  322.         zcopy(z1, res);
  323.         return;
  324.     }
  325.  
  326.     /*
  327.      * Now make sure the input number is within the modulo range.
  328.      * If not, then reduce it to be within range and make the
  329.      * quick check again.
  330.      */
  331.     sign = z1.sign;
  332.     z1.sign = 0;
  333.     cv = zrel(z1, z2);
  334.     if (cv == 0) {
  335.         *res = _zero_;
  336.         return;
  337.     }
  338.     tmp1 = z1;
  339.     if (cv > 0) {
  340.         z1.sign = (BOOL)sign;
  341.         zmod(z1, z2, &tmp1);
  342.         if (tmp1.len < z2.len - 1) {
  343.             *res = tmp1;
  344.             return;
  345.         }
  346.         sign = 0;
  347.     }
  348.  
  349.     /*
  350.      * Now calculate the difference of the modulus and the absolute
  351.      * value of the original number.  Compare the original number with
  352.      * the difference, and return the one with the smallest absolute
  353.      * value, with the correct sign.  If the two values are equal, then
  354.      * return the positive result.
  355.      */
  356.     zsub(z2, tmp1, &tmp2);
  357.     cv = zrel(tmp1, tmp2);
  358.     if (cv < 0) {
  359.         freeh(tmp2.v);
  360.         tmp1.sign = (BOOL)sign;
  361.         if (tmp1.v == z1.v)
  362.             zcopy(tmp1, res);
  363.         else
  364.             *res = tmp1;
  365.     } else {
  366.         if (cv)
  367.             tmp2.sign = !sign;
  368.         if (tmp1.v != z1.v)
  369.             freeh(tmp1.v);
  370.         *res = tmp2;
  371.     }
  372. }
  373.  
  374.  
  375. /*
  376.  * Compare two numbers for equality modulo a third number.
  377.  * The two numbers to be compared can be negative or out of modulo range.
  378.  * Returns TRUE if the numbers are not congruent, and FALSE if they are
  379.  * congruent.
  380.  */
  381. BOOL
  382. zcmpmod(z1, z2, z3)
  383.     ZVALUE z1;        /* first number to be compared */
  384.     ZVALUE z2;        /* second number to be compared */
  385.     ZVALUE z3;        /* modulus */
  386. {
  387.     ZVALUE tmp1, tmp2, tmp3;
  388.     FULL digit;
  389.     LEN len;
  390.     int cv;
  391.  
  392.     if (isneg(z3) || iszero(z3))
  393.         error("Non-positive modulus in zcmpmod");
  394.     if (istwo(z3))
  395.         return (((z1.v[0] + z2.v[0]) & 0x1) != 0);
  396.  
  397.     /*
  398.      * If the two numbers are equal, then their mods are equal.
  399.      */
  400.     if ((z1.sign == z2.sign) && (z1.len == z2.len) &&
  401.         (z1.v[0] == z2.v[0]) && (zcmp(z1, z2) == 0))
  402.             return FALSE;
  403.  
  404.     /*
  405.      * If both numbers are negative, then we can make them positive.
  406.      */
  407.     if (isneg(z1) && isneg(z2)) {
  408.         z1.sign = 0;
  409.         z2.sign = 0;
  410.     }
  411.  
  412.     /*
  413.      * For small negative numbers, make them positive before comparing.
  414.      * In any case, the resulting numbers are in tmp1 and tmp2.
  415.      */
  416.     tmp1 = z1;
  417.     tmp2 = z2;
  418.     len = z3.len;
  419.     digit = z3.v[len - 1];
  420.  
  421.     if (isneg(z1) && ((z1.len < len) ||
  422.         ((z1.len == len) && (z1.v[z1.len - 1] < digit))))
  423.             zadd(z1, z3, &tmp1);
  424.  
  425.     if (isneg(z2) && ((z2.len < len) ||
  426.         ((z2.len == len) && (z2.v[z2.len - 1] < digit))))
  427.             zadd(z2, z3, &tmp2);
  428.  
  429.     /*
  430.      * Now compare the two numbers for equality.
  431.      * If they are equal we are all done.
  432.      */
  433.     if (zcmp(tmp1, tmp2) == 0) {
  434.         if (tmp1.v != z1.v)
  435.             freeh(tmp1.v);
  436.         if (tmp2.v != z2.v)
  437.             freeh(tmp2.v);
  438.         return FALSE;
  439.     }
  440.  
  441.     /*
  442.      * They are not identical.  Now if both numbers are positive
  443.      * and less than the modulus, then they are definitely not equal.
  444.      */
  445.     if ((tmp1.sign == tmp2.sign) &&
  446.         ((tmp1.len < len) || (zrel(tmp1, z3) < 0)) &&
  447.         ((tmp2.len < len) || (zrel(tmp2, z3) < 0)))
  448.     {
  449.         if (tmp1.v != z1.v)
  450.             freeh(tmp1.v);
  451.         if (tmp2.v != z2.v)
  452.             freeh(tmp2.v);
  453.         return TRUE;
  454.     }
  455.  
  456.     /*
  457.      * Either one of the numbers is negative or is large.
  458.      * So do the standard thing and subtract the two numbers.
  459.      * Then they are equal if the result is 0 (mod z3).
  460.      */
  461.     zsub(tmp1, tmp2, &tmp3);
  462.     if (tmp1.v != z1.v)
  463.         freeh(tmp1.v);
  464.     if (tmp2.v != z2.v)
  465.         freeh(tmp2.v);
  466.  
  467.     /*
  468.      * Compare the result with the modulus to see if it is equal to
  469.      * or less than the modulus.  If so, we know the mod result.
  470.      */
  471.     tmp3.sign = 0;
  472.     cv = zrel(tmp3, z3);
  473.     if (cv == 0) {
  474.         freeh(tmp3.v);
  475.         return FALSE;
  476.     }
  477.     if (cv < 0) {
  478.         freeh(tmp3.v);
  479.         return TRUE;
  480.     }
  481.  
  482.     /*
  483.      * We are forced to actually do the division.
  484.      * The numbers are congruent if the result is zero.
  485.      */
  486.     zmod(tmp3, z3, &tmp1);
  487.     freeh(tmp3.v);
  488.     if (iszero(tmp1)) {
  489.         freeh(tmp1.v);
  490.         return FALSE;
  491.     } else {
  492.         freeh(tmp1.v);
  493.         return TRUE;
  494.     }
  495. }
  496.  
  497.  
  498. /*
  499.  * Compute the result of raising one number to a power modulo another number.
  500.  * That is, this computes:  a^b (modulo c).
  501.  * This calculates the result by examining the power POWBITS bits at a time,
  502.  * using a small table of POWNUMS low powers to calculate powers for those bits,
  503.  * and repeated squaring and multiplying by the partial powers to generate
  504.  * the complete power.  If the power being raised to is high enough, then
  505.  * this uses the REDC algorithm to avoid doing many divisions.  When using
  506.  * REDC, multiple calls to this routine using the same modulus will be
  507.  * slightly faster.
  508.  */
  509. void
  510. zpowermod(z1, z2, z3, res)
  511.     ZVALUE z1, z2, z3, *res;
  512. {
  513.     HALF *hp;        /* pointer to current word of the power */
  514.     REDC *rp;        /* REDC information to be used */
  515.     ZVALUE *pp;        /* pointer to low power table */
  516.     ZVALUE ans, temp;    /* calculation values */
  517.     ZVALUE modpow;        /* current small power */
  518.     ZVALUE lowpowers[POWNUMS];    /* low powers */
  519.     int sign;        /* original sign of number */
  520.     int curshift;        /* shift value for word of power */
  521.     HALF curhalf;        /* current word of power */
  522.     unsigned int curpow;    /* current low power */
  523.     unsigned int curbit;    /* current bit of low power */
  524.     int i;
  525.  
  526.     if (isneg(z3) || iszero(z3))
  527.         error("Non-positive modulus in zpowermod");
  528.     if (isneg(z2))
  529.         error("Negative power in zpowermod");
  530.  
  531.     sign = z1.sign;
  532.     z1.sign = 0;
  533.  
  534.     /*
  535.      * Check easy cases first.
  536.      */
  537.     if (iszero(z1) || isunit(z3)) {        /* 0^x or mod 1 */
  538.         *res = _zero_;
  539.         return;
  540.     }
  541.     if (istwo(z3)) {            /* mod 2 */
  542.         if (isodd(z1))
  543.             *res = _one_;
  544.         else
  545.             *res = _zero_;
  546.         return;
  547.     }
  548.     if (isunit(z1) && (!sign || iseven(z2))) {
  549.         /* 1^x or (-1)^(2x) */
  550.         *res = _one_;
  551.         return;
  552.     }
  553.  
  554.     /*
  555.      * Normalize the number being raised to be non-negative and to lie
  556.      * within the modulo range.  Then check for zero or one specially.
  557.      */
  558.     zmod(z1, z3, &temp);
  559.     if (iszero(temp)) {
  560.         freeh(temp.v);
  561.         *res = _zero_;
  562.         return;
  563.     }
  564.     z1 = temp;
  565.     if (sign) {
  566.         zsub(z3, z1, &temp);
  567.         freeh(z1.v);
  568.         z1 = temp;
  569.     }
  570.     if (isunit(z1)) {
  571.         freeh(z1.v);
  572.         *res = _one_;
  573.         return;
  574.     }
  575.  
  576.     /*
  577.      * If the modulus is odd, large enough, is not one less than an
  578.      * exact power of two, and if the power is large enough, then use
  579.      * the REDC algorithm.  The size where this is done is configurable.
  580.      */
  581.     if ((z2.len > 1) && (z3.len >= _pow2_) && isodd(z3)
  582.         && !zisallbits(z3))
  583.     {
  584.         if (powermodredc && zcmp(powermodredc->mod, z3)) {
  585.             zredcfree(powermodredc);
  586.             powermodredc = NULL;
  587.         }
  588.         if (powermodredc == NULL)
  589.             powermodredc = zredcalloc(z3);
  590.         rp = powermodredc;
  591.         zredcencode(rp, z1, &temp);
  592.         zredcpower(rp, temp, z2, &z1);
  593.         freeh(temp.v);
  594.         zredcdecode(rp, z1, res);
  595.         freeh(z1.v);
  596.         return;
  597.     }
  598.  
  599.     /*
  600.      * Modulus or power is small enough to perform the power raising
  601.      * directly.  Initialize the table of powers.
  602.      */
  603.     for (pp = &lowpowers[2]; pp < &lowpowers[POWNUMS]; pp++)
  604.         pp->len = 0;
  605.     lowpowers[0] = _one_;
  606.     lowpowers[1] = z1;
  607.     ans = _one_;
  608.  
  609.     hp = &z2.v[z2.len - 1];
  610.     curhalf = *hp;
  611.     curshift = BASEB - POWBITS;
  612.     while (curshift && ((curhalf >> curshift) == 0))
  613.         curshift -= POWBITS;
  614.  
  615.     /*
  616.      * Calculate the result by examining the power POWBITS bits at a time,
  617.      * and use the table of low powers at each iteration.
  618.      */
  619.     for (;;) {
  620.         curpow = (curhalf >> curshift) & (POWNUMS - 1);
  621.         pp = &lowpowers[curpow];
  622.  
  623.         /*
  624.          * If the small power is not yet saved in the table, then
  625.          * calculate it and remember it in the table for future use.
  626.          */
  627.         if (pp->len == 0) {
  628.             if (curpow & 0x1)
  629.                 zcopy(z1, &modpow);
  630.             else
  631.                 modpow = _one_;
  632.  
  633.             for (curbit = 0x2; curbit <= curpow; curbit *= 2) {
  634.                 pp = &lowpowers[curbit];
  635.                 if (pp->len == 0) {
  636.                     zsquare(lowpowers[curbit/2], &temp);
  637.                     zmod(temp, z3, pp);
  638.                     freeh(temp.v);
  639.                 }
  640.                 if (curbit & curpow) {
  641.                     zmul(*pp, modpow, &temp);
  642.                     freeh(modpow.v);
  643.                     zmod(temp, z3, &modpow);
  644.                     freeh(temp.v);
  645.                 }
  646.             }
  647.             pp = &lowpowers[curpow];
  648.             *pp = modpow;
  649.         }
  650.  
  651.         /*
  652.          * If the power is nonzero, then accumulate the small power
  653.          * into the result.
  654.          */
  655.         if (curpow) {
  656.             zmul(ans, *pp, &temp);
  657.             freeh(ans.v);
  658.             zmod(temp, z3, &ans);
  659.             freeh(temp.v);
  660.         }
  661.  
  662.         /*
  663.          * Select the next POWBITS bits of the power, if there is
  664.          * any more to generate.
  665.          */
  666.         curshift -= POWBITS;
  667.         if (curshift < 0) {
  668.             if (hp-- == z2.v)
  669.                 break;
  670.             curhalf = *hp;
  671.             curshift = BASEB - POWBITS;
  672.         }
  673.  
  674.         /*
  675.          * Square the result POWBITS times to make room for the next
  676.          * chunk of bits.
  677.          */
  678.         for (i = 0; i < POWBITS; i++) {
  679.             zsquare(ans, &temp);
  680.             freeh(ans.v);
  681.             zmod(temp, z3, &ans);
  682.             freeh(temp.v);
  683.         }
  684.     }
  685.  
  686.     for (pp = &lowpowers[2]; pp < &lowpowers[POWNUMS]; pp++) {
  687.         if (pp->len)
  688.             freeh(pp->v);
  689.     }
  690.     *res = ans;
  691. }
  692.  
  693.  
  694. /*
  695.  * Initialize the REDC algorithm for a particular modulus,
  696.  * returning a pointer to a structure that is used for other
  697.  * REDC calls.  An error is generated if the structure cannot
  698.  * be allocated.  The modulus must be odd and positive.
  699.  */
  700. REDC *
  701. zredcalloc(z1)
  702.     ZVALUE z1;        /* modulus to initialize for */
  703. {
  704.     REDC *rp;        /* REDC information */
  705.     ZVALUE tmp;
  706.     long bit;
  707.  
  708.     if (iseven(z1) || isneg(z1))
  709.         error("REDC requires positive odd modulus");
  710.  
  711.     rp = (REDC *) malloc(sizeof(REDC));
  712.     if (rp == NULL)
  713.         error("Cannot allocate REDC structure");
  714.  
  715.     /*
  716.      * Round up the binary modulus to the next power of two
  717.      * which is at a word boundary.  Then the shift and modulo
  718.      * operations mod the binary modulus can be done very cheaply.
  719.      * Calculate the REDC format for the number 1 for future use.
  720.      */
  721.     bit = zhighbit(z1) + 1;
  722.     if (bit % BASEB)
  723.         bit += (BASEB - (bit % BASEB));
  724.     zcopy(z1, &rp->mod);
  725.     zbitvalue(bit, &tmp);
  726.     z1.sign = 1;
  727.     (void) zmodinv(z1, tmp, &rp->inv);
  728.     zmod(tmp, rp->mod, &rp->one);
  729.     freeh(tmp.v);
  730.     rp->len = bit / BASEB;
  731.     return rp;
  732. }
  733.  
  734.  
  735. /*
  736.  * Free any numbers associated with the specified REDC structure,
  737.  * and then the REDC structure itself.
  738.  */
  739. void
  740. zredcfree(rp)
  741.     REDC *rp;        /* REDC information to be cleared */
  742. {
  743.     freeh(rp->mod.v);
  744.     freeh(rp->inv.v);
  745.     freeh(rp->one.v);
  746.     free(rp);
  747. }
  748.  
  749.  
  750. /*
  751.  * Convert a normal number into the specified REDC format.
  752.  * The number to be converted can be negative or out of modulo range.
  753.  * The resulting number can be used for multiplying, adding, subtracting,
  754.  * or comparing with any other such converted numbers, as if the numbers
  755.  * were being calculated modulo the number which initialized the REDC
  756.  * information.  When the final value is unconverted, the result is the
  757.  * same as if the usual operations were done with the original numbers.
  758.  */
  759. void
  760. zredcencode(rp, z1, res)
  761.     REDC *rp;        /* REDC information */
  762.     ZVALUE z1;        /* number to be converted */
  763.     ZVALUE *res;        /* returned converted number */
  764. {
  765.     ZVALUE tmp1, tmp2;
  766.  
  767.     /*
  768.      * Handle the cases 0, 1, -1, and 2 specially since these are
  769.      * easy to calculate.  Zero transforms to zero, and the others
  770.      * can be obtained from the precomputed REDC format for 1 since
  771.      * addition and subtraction act normally for REDC format numbers.
  772.      */
  773.     if (iszero(z1)) {
  774.         *res = _zero_;
  775.         return;
  776.     }
  777.     if (isone(z1)) {
  778.         zcopy(rp->one, res);
  779.         return;
  780.     }
  781.     if (isunit(z1)) {
  782.         zsub(rp->mod, rp->one, res);
  783.         return;
  784.     }
  785.     if (istwo(z1)) {
  786.         zadd(rp->one, rp->one, &tmp1);
  787.         if (zrel(tmp1, rp->mod) < 0) {
  788.             *res = tmp1;
  789.             return;
  790.         }
  791.         zsub(tmp1, rp->mod, res);
  792.         freeh(tmp1.v);
  793.         return;
  794.     }
  795.  
  796.     /*
  797.      * Not a trivial number to convert, so do the full transformation.
  798.      * Convert negative numbers to positive numbers before converting.
  799.      */
  800.     tmp1.len = 0;
  801.     if (isneg(z1)) {
  802.         zmod(z1, rp->mod, &tmp1);
  803.         z1 = tmp1;
  804.     }
  805.     zshift(z1, rp->len * BASEB, &tmp2);
  806.     if (tmp1.len)
  807.         freeh(tmp1.v);
  808.     zmod(tmp2, rp->mod, res);
  809.     freeh(tmp2.v);
  810. }
  811.  
  812.  
  813. /*
  814.  * The REDC algorithm used to convert numbers out of REDC format and also
  815.  * used after multiplication of two REDC numbers.  Using this routine
  816.  * avoids any divides, replacing the divide by two multiplications.
  817.  * If the numbers are very large, then these two multiplies will be
  818.  * quicker than the divide, since dividing is harder than multiplying.
  819.  */
  820. void
  821. zredcdecode(rp, z1, res)
  822.     REDC *rp;        /* REDC information */
  823.     ZVALUE z1;        /* number to be transformed */
  824.     ZVALUE *res;        /* returned transformed number */
  825. {
  826.     ZVALUE tmp1, tmp2;    /* temporaries */
  827.     HALF *hp;        /* saved pointer to tmp2 value */
  828.  
  829.     if (isneg(z1))
  830.         error("Negative number for zredc");
  831.  
  832.     /*
  833.      * Check first for the special values for 0 and 1 that are easy.
  834.      */
  835.     if (iszero(z1)) {
  836.         *res = _zero_;
  837.         return;
  838.     }
  839.     if ((z1.len == rp->one.len) && (z1.v[0] == rp->one.v[0]) &&
  840.         (zcmp(z1, rp->one) == 0)) {
  841.             *res = _one_;
  842.             return;
  843.     }
  844.  
  845.     /*
  846.      * First calculate the following:
  847.      *     tmp2 = ((z1 % 2^bitnum) * inv) % 2^bitnum.
  848.      * The mod operations can be done with no work since the bit
  849.      * number was selected as a multiple of the word size.  Just
  850.      * reduce the sizes of the numbers as required.
  851.      */
  852.     tmp1 = z1;
  853.     if (tmp1.len > rp->len)
  854.         tmp1.len = rp->len;
  855.     zmul(tmp1, rp->inv, &tmp2);
  856.     if (tmp2.len > rp->len)
  857.         tmp2.len = rp->len;
  858.  
  859.     /*
  860.      * Next calculate the following:
  861.      *    res = (z1 + tmp2 * modulus) / 2^bitnum
  862.      * The division by a power of 2 is always exact, and requires no
  863.      * work.  Just adjust the address and length of the number to do
  864.      * the divide, but save the original pointer for freeing later.
  865.      */
  866.     zmul(tmp2, rp->mod, &tmp1);
  867.     freeh(tmp2.v);
  868.     zadd(z1, tmp1, &tmp2);
  869.     freeh(tmp1.v);
  870.     hp = tmp2.v;
  871.     if (tmp2.len <= rp->len) {
  872.         freeh(hp);
  873.         *res = _zero_;
  874.         return;
  875.     }
  876.     tmp2.v += rp->len;
  877.     tmp2.len -= rp->len;
  878.  
  879.     /*
  880.      * Finally do a final modulo by a simple subtraction if necessary.
  881.      * This is all that is needed because the previous calculation is
  882.      * guaranteed to always be less than twice the modulus.
  883.      */
  884.     if (zrel(tmp2, rp->mod) < 0)
  885.         zcopy(tmp2, res);
  886.     else
  887.         zsub(tmp2, rp->mod, res);
  888.     freeh(hp);
  889. }
  890.  
  891.  
  892. /*
  893.  * Multiply two numbers in REDC format together producing a result also
  894.  * in REDC format.  If the result is converted back to a normal number,
  895.  * then the result is the same as the modulo'd multiplication of the
  896.  * original numbers before they were converted to REDC format.  This
  897.  * calculation is done in one of two ways, depending on the size of the
  898.  * modulus.  For large numbers, the REDC definition is used directly
  899.  * which involves three multiplies overall.  For small numbers, a
  900.  * complicated routine is used which does the indicated multiplication
  901.  * and the REDC algorithm at the same time to produce the result.
  902.  */
  903. void
  904. zredcmul(rp, z1, z2, res)
  905.     REDC *rp;        /* REDC information */
  906.     ZVALUE z1;        /* first REDC number to be multiplied */
  907.     ZVALUE z2;        /* second REDC number to be multiplied */
  908.     ZVALUE *res;        /* resulting REDC number */
  909. {
  910.     FULL mulb;
  911.     FULL muln;
  912.     HALF *h1;
  913.     HALF *h2;
  914.     HALF *h3;
  915.     HALF *hd;
  916.     HALF Ninv;
  917.     HALF topdigit;
  918.     LEN modlen;
  919.     LEN len;
  920.     LEN len2;
  921.     SIUNION sival1;
  922.     SIUNION sival2;
  923.     SIUNION sival3;
  924.     SIUNION carry;
  925.     ZVALUE tmp;
  926.  
  927.     if (isneg(z1) || (z1.len > rp->mod.len) ||
  928.         isneg(z2) || (z2.len > rp->mod.len))
  929.             error("Negative or too large number in zredcmul");
  930.  
  931.     /*
  932.      * Check for special values which we easily know the answer.
  933.      */
  934.     if (iszero(z1) || iszero(z2)) {
  935.         *res = _zero_;
  936.         return;
  937.     }
  938.  
  939.     if ((z1.len == rp->one.len) && (z1.v[0] == rp->one.v[0]) &&
  940.         (zcmp(z1, rp->one) == 0)) {
  941.             zcopy(z2, res);
  942.             return;
  943.     }
  944.  
  945.     if ((z2.len == rp->one.len) && (z2.v[0] == rp->one.v[0]) &&
  946.         (zcmp(z2, rp->one) == 0)) {
  947.             zcopy(z1, res);
  948.             return;
  949.     }
  950.  
  951.     /*
  952.      * If the size of the modulus is large, then just do the multiply,
  953.      * followed by the two multiplies contained in the REDC routine.
  954.      * This will be quicker than directly doing the REDC calculation
  955.      * because of the O(N^1.585) speed of the multiplies.  The size
  956.      * of the number which this is done is configurable.
  957.      */
  958.     if (rp->mod.len >= _redc2_) {
  959.         zmul(z1, z2, &tmp);
  960.         zredcdecode(rp, tmp, res);
  961.         freeh(tmp.v);
  962.         return;
  963.     }
  964.  
  965.     /*
  966.      * The number is small enough to calculate by doing the O(N^2) REDC
  967.      * algorithm directly.  This algorithm performs the multiplication and
  968.      * the reduction at the same time.  Notice the obscure facts that
  969.      * only the lowest word of the inverse value is used, and that
  970.      * there is no shifting of the partial products as there is in a
  971.      * normal multiply.
  972.      */
  973.     modlen = rp->mod.len;
  974.     Ninv = rp->inv.v[0];
  975.  
  976.     /*
  977.      * Allocate the result and clear it.
  978.      * The size of the result will be equal to or smaller than
  979.      * the modulus size.
  980.      */
  981.     res->sign = 0;
  982.     res->len = modlen;
  983.     res->v = alloc(modlen);
  984.  
  985.     hd = res->v;
  986.     len = modlen;
  987.     while (len--)
  988.         *hd++ = 0;
  989.  
  990.     /*
  991.      * Do this outermost loop over all the digits of z1.
  992.      */
  993.     h1 = z1.v;
  994.     len = z1.len;
  995.     while (len--) {
  996.         /*
  997.          * Start off with the next digit of z1, the first
  998.          * digit of z2, and the first digit of the modulus.
  999.          */
  1000.         mulb = (FULL) *h1++;
  1001.         h2 = z2.v;
  1002.         h3 = rp->mod.v;
  1003.         hd = res->v;
  1004.         sival1.ivalue = mulb * ((FULL) *h2++) + ((FULL) *hd++);
  1005.         muln = ((HALF) (sival1.silow * Ninv));
  1006.         sival2.ivalue = muln * ((FULL) *h3++);
  1007.         sival3.ivalue = ((FULL) sival1.silow) + ((FULL) sival2.silow);
  1008.         carry.ivalue = ((FULL) sival1.sihigh) + ((FULL) sival2.sihigh)
  1009.             + ((FULL) sival3.sihigh);
  1010.  
  1011.         /*
  1012.          * Do this innermost loop for each digit of z2, except
  1013.          * for the first digit which was just done above.
  1014.          */
  1015.         len2 = z2.len;
  1016.         while (--len2 > 0) {
  1017.             sival1.ivalue = mulb * ((FULL) *h2++);
  1018.             sival2.ivalue = muln * ((FULL) *h3++);
  1019.             sival3.ivalue = ((FULL) sival1.silow)
  1020.                 + ((FULL) sival2.silow)
  1021.                 + ((FULL) *hd)
  1022.                 + ((FULL) carry.silow);
  1023.             carry.ivalue = ((FULL) sival1.sihigh)
  1024.                 + ((FULL) sival2.sihigh)
  1025.                 + ((FULL) sival3.sihigh)
  1026.                 + ((FULL) carry.sihigh);
  1027.  
  1028.             hd[-1] = sival3.silow;
  1029.             hd++;
  1030.         }
  1031.  
  1032.         /*
  1033.          * Now continue the loop as necessary so the total number
  1034.          * of interations is equal to the size of the modulus.
  1035.          * This acts as if the innermost loop was repeated for
  1036.          * high digits of z2 that are zero.
  1037.          */
  1038.         len2 = modlen - z2.len;
  1039.         while (len2--) {
  1040.             sival2.ivalue = muln * ((FULL) *h3++);
  1041.             sival3.ivalue = ((FULL) sival2.silow)
  1042.                 + ((FULL) *hd)
  1043.                 + ((FULL) carry.silow);
  1044.             carry.ivalue = ((FULL) sival2.sihigh)
  1045.                 + ((FULL) sival3.sihigh)
  1046.                 + ((FULL) carry.sihigh);
  1047.  
  1048.             hd[-1] = sival3.silow;
  1049.             hd++;
  1050.         }
  1051.  
  1052.         res->v[modlen - 1] = carry.silow;
  1053.         topdigit = carry.sihigh;
  1054.     }
  1055.  
  1056.     /*
  1057.      * Now continue the loop as necessary so the total number
  1058.      * of interations is equal to the size of the modulus.
  1059.      * This acts as if the outermost loop was repeated for high
  1060.      * digits of z1 that are zero.
  1061.      */
  1062.     len = modlen - z1.len;
  1063.     while (len--) {
  1064.         /*
  1065.          * Start off with the first digit of the modulus.
  1066.          */
  1067.         h3 = rp->mod.v;
  1068.         hd = res->v;
  1069.         muln = ((HALF) (*hd * Ninv));
  1070.         sival2.ivalue = muln * ((FULL) *h3++);
  1071.         sival3.ivalue = ((FULL) *hd++) + ((FULL) sival2.silow);
  1072.         carry.ivalue = ((FULL) sival2.sihigh) + ((FULL) sival3.sihigh);
  1073.  
  1074.         /*
  1075.          * Do this innermost loop for each digit of the modulus,
  1076.          * except for the first digit which was just done above.
  1077.          */
  1078.         len2 = modlen;
  1079.         while (--len2 > 0) {
  1080.             sival2.ivalue = muln * ((FULL) *h3++);
  1081.             sival3.ivalue = ((FULL) sival2.silow)
  1082.                 + ((FULL) *hd)
  1083.                 + ((FULL) carry.silow);
  1084.             carry.ivalue = ((FULL) sival2.sihigh)
  1085.                 + ((FULL) sival3.sihigh)
  1086.                 + ((FULL) carry.sihigh);
  1087.  
  1088.             hd[-1] = sival3.silow;
  1089.             hd++;
  1090.         }
  1091.         res->v[modlen - 1] = carry.silow;
  1092.         topdigit = carry.sihigh;
  1093.     }
  1094.  
  1095.     /*
  1096.      * Determine the true size of the result, taking the top digit of
  1097.      * the current result into account.  The top digit is not stored in
  1098.      * the number because it is temporary and would become zero anyway
  1099.      * after the final subtraction is done.
  1100.      */
  1101.     if (topdigit == 0) {
  1102.         len = modlen;
  1103.         hd = &res->v[len - 1];
  1104.         while ((*hd == 0) && (len > 1)) {
  1105.             hd--;
  1106.             len--;
  1107.         }
  1108.         res->len = len;
  1109.     }
  1110.  
  1111.     /*
  1112.      * Compare the result with the modulus.
  1113.      * If it is less than the modulus, then the calculation is complete.
  1114.      */
  1115.     if ((topdigit == 0) && ((len < modlen)
  1116.         || (res->v[len - 1] < rp->mod.v[len - 1])
  1117.         || (zrel(*res, rp->mod) < 0)))
  1118.             return;
  1119.  
  1120.     /*
  1121.      * Do a subtraction to reduce the result to a value less than
  1122.      * the modulus.  The REDC algorithm guarantees that a single subtract
  1123.      * is all that is needed.  Ignore any borrowing from the possible
  1124.      * highest word of the current result because that would affect
  1125.      * only the top digit value that was not stored and would become
  1126.      * zero anyway.
  1127.      */
  1128.     carry.ivalue = 0;
  1129.     h1 = rp->mod.v;
  1130.     hd = res->v;
  1131.     len = modlen;
  1132.     while (len--) {
  1133.         carry.ivalue = BASE1 - ((FULL) *hd) + ((FULL) *h1++)
  1134.             + ((FULL) carry.silow);
  1135.         *hd++ = BASE1 - carry.silow;
  1136.         carry.silow = carry.sihigh;
  1137.     }
  1138.  
  1139.     /*
  1140.      * Now finally recompute the size of the result.
  1141.      */
  1142.     len = modlen;
  1143.     hd = &res->v[len - 1];
  1144.     while ((*hd == 0) && (len > 1)) {
  1145.         hd--;
  1146.         len--;
  1147.     }
  1148.     res->len = len;
  1149. }
  1150.  
  1151.  
  1152. /*
  1153.  * Square a number in REDC format producing a result also in REDC format.
  1154.  */
  1155. void
  1156. zredcsquare(rp, z1, res)
  1157.     REDC *rp;        /* REDC information */
  1158.     ZVALUE z1;        /* REDC number to be squared */
  1159.     ZVALUE *res;        /* resulting REDC number */
  1160. {
  1161.     ZVALUE tmp;
  1162.  
  1163.     if (isneg(z1))
  1164.         error("Negative number in zredcsquare");
  1165.     if (iszero(z1)) {
  1166.         *res = _zero_;
  1167.         return;
  1168.     }
  1169.     if ((z1.len == rp->one.len) && (z1.v[0] == rp->one.v[0]) &&
  1170.         (zcmp(z1, rp->one) == 0)) {
  1171.             zcopy(z1, res);
  1172.             return;
  1173.     }
  1174.  
  1175.     /*
  1176.      * If the modulus is small enough, then call the multiply
  1177.      * routine to produce the result.  Otherwise call the O(N^1.585)
  1178.      * routines to get the answer.
  1179.      */
  1180.     if (rp->mod.len < _redc2_) {
  1181.         zredcmul(rp, z1, z1, res);
  1182.         return;
  1183.     }
  1184.     zsquare(z1, &tmp);
  1185.     zredcdecode(rp, tmp, res);
  1186.     freeh(tmp.v);
  1187. }
  1188.  
  1189.  
  1190. /*
  1191.  * Compute the result of raising a REDC format number to a power.
  1192.  * The result is within the range 0 to the modulus - 1.
  1193.  * This calculates the result by examining the power POWBITS bits at a time,
  1194.  * using a small table of POWNUMS low powers to calculate powers for those bits,
  1195.  * and repeated squaring and multiplying by the partial powers to generate
  1196.  * the complete power.
  1197.  */
  1198. void
  1199. zredcpower(rp, z1, z2, res)
  1200.     REDC *rp;        /* REDC information */
  1201.     ZVALUE z1;        /* REDC number to be raised */
  1202.     ZVALUE z2;        /* normal number to raise number to */
  1203.     ZVALUE *res;        /* result */
  1204. {
  1205.     HALF *hp;        /* pointer to current word of the power */
  1206.     ZVALUE *pp;        /* pointer to low power table */
  1207.     ZVALUE ans, temp;    /* calculation values */
  1208.     ZVALUE modpow;        /* current small power */
  1209.     ZVALUE lowpowers[POWNUMS];    /* low powers */
  1210.     int curshift;        /* shift value for word of power */
  1211.     HALF curhalf;        /* current word of power */
  1212.     unsigned int curpow;    /* current low power */
  1213.     unsigned int curbit;    /* current bit of low power */
  1214.     int i;
  1215.  
  1216.     if (isneg(z1))
  1217.         error("Negative number in zredcpower");
  1218.     if (isneg(z2))
  1219.         error("Negative power in zredcpower");
  1220.  
  1221.     /*
  1222.      * Check for zero or the REDC format for one.
  1223.      */
  1224.     if (iszero(z1) || isunit(rp->mod)) {
  1225.         *res = _zero_;
  1226.         return;
  1227.     }
  1228.     if (zcmp(z1, rp->one) == 0) {
  1229.         zcopy(rp->one, res);
  1230.         return;
  1231.     }
  1232.  
  1233.     /*
  1234.      * See if the number being raised is the REDC format for -1.
  1235.      * If so, then the answer is the REDC format for one or minus one.
  1236.      * To do this check, calculate the REDC format for -1.
  1237.      */
  1238.     if (((HALF)(z1.v[0] + rp->one.v[0])) == rp->mod.v[0]) {
  1239.         zsub(rp->mod, rp->one, &temp);
  1240.         if (zcmp(z1, temp) == 0) {
  1241.             if (isodd(z2)) {
  1242.                 *res = temp;
  1243.                 return;
  1244.             }
  1245.             freeh(temp.v);
  1246.             zcopy(rp->one, res);
  1247.             return;
  1248.         }
  1249.         freeh(temp.v);
  1250.     }
  1251.  
  1252.     for (pp = &lowpowers[2]; pp < &lowpowers[POWNUMS]; pp++)
  1253.         pp->len = 0;
  1254.     zcopy(rp->one, &lowpowers[0]);
  1255.     zcopy(z1, &lowpowers[1]);
  1256.     zcopy(rp->one, &ans);
  1257.  
  1258.     hp = &z2.v[z2.len - 1];
  1259.     curhalf = *hp;
  1260.     curshift = BASEB - POWBITS;
  1261.     while (curshift && ((curhalf >> curshift) == 0))
  1262.         curshift -= POWBITS;
  1263.  
  1264.     /*
  1265.      * Calculate the result by examining the power POWBITS bits at a time,
  1266.      * and use the table of low powers at each iteration.
  1267.      */
  1268.     for (;;) {
  1269.         curpow = (curhalf >> curshift) & (POWNUMS - 1);
  1270.         pp = &lowpowers[curpow];
  1271.  
  1272.         /*
  1273.          * If the small power is not yet saved in the table, then
  1274.          * calculate it and remember it in the table for future use.
  1275.          */
  1276.         if (pp->len == 0) {
  1277.             if (curpow & 0x1)
  1278.                 zcopy(z1, &modpow);
  1279.             else
  1280.                 zcopy(rp->one, &modpow);
  1281.  
  1282.             for (curbit = 0x2; curbit <= curpow; curbit *= 2) {
  1283.                 pp = &lowpowers[curbit];
  1284.                 if (pp->len == 0)
  1285.                     zredcsquare(rp, lowpowers[curbit/2],
  1286.                         pp);
  1287.                 if (curbit & curpow) {
  1288.                     zredcmul(rp, *pp, modpow, &temp);
  1289.                     freeh(modpow.v);
  1290.                     modpow = temp;
  1291.                 }
  1292.             }
  1293.             pp = &lowpowers[curpow];
  1294.             *pp = modpow;
  1295.         }
  1296.  
  1297.         /*
  1298.          * If the power is nonzero, then accumulate the small power
  1299.          * into the result.
  1300.          */
  1301.         if (curpow) {
  1302.             zredcmul(rp, ans, *pp, &temp);
  1303.             freeh(ans.v);
  1304.             ans = temp;
  1305.         }
  1306.  
  1307.         /*
  1308.          * Select the next POWBITS bits of the power, if there is
  1309.          * any more to generate.
  1310.          */
  1311.         curshift -= POWBITS;
  1312.         if (curshift < 0) {
  1313.             if (hp-- == z2.v)
  1314.                 break;
  1315.             curhalf = *hp;
  1316.             curshift = BASEB - POWBITS;
  1317.         }
  1318.  
  1319.         /*
  1320.          * Square the result POWBITS times to make room for the next
  1321.          * chunk of bits.
  1322.          */
  1323.         for (i = 0; i < POWBITS; i++) {
  1324.             zredcsquare(rp, ans, &temp);
  1325.             freeh(ans.v);
  1326.             ans = temp;
  1327.         }
  1328.     }
  1329.  
  1330.     for (pp = lowpowers; pp < &lowpowers[POWNUMS]; pp++) {
  1331.         if (pp->len)
  1332.             freeh(pp->v);
  1333.     }
  1334.     *res = ans;
  1335. }
  1336.  
  1337. /* END CODE */
  1338.